In [4]:
import sys

# Our numerical workhorses
import numpy as np
import pandas as pd
import scipy.optimize
from scipy.integrate import quad
import scipy.stats as st

# Numerical differentiation package
import numdifftools as ndt

# MCMC package
import emcee

# MCMC results visualization package
import corner

# Import pyplot for plotting
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import LightSource

# Seaborn, useful for graphics
#import seaborn as sns

# Bokeh stuff
import bokeh.charts
import bokeh.io

# Magic function to make matplotlib inline; other style specs must come AFTER
%matplotlib inline

# This enables high res graphics inline (only use with static plots (non-Bokeh))
# SVG is preferred, but there is a bug in Jupyter with vertical lines
%config InlineBackend.figure_formats = {'svg', 'retina'}

# JB's favorite Seaborn settings for notebooks
rc = {'lines.linewidth': 2, 
      'axes.labelsize': 18, 
      'axes.titlesize': 18, 
      'axes.facecolor': 'DFDFE5'}
#sns.set_context('notebook', font_scale=1.5, rc={'lines.linewidth': 2.5})
#sns.set_style('darkgrid', {'axes.facecolor': '(0.875, 0.875, 0.9)'})

# Display graphics in this notebook
bokeh.io.output_notebook()
Loading BokehJS ...
In [5]:
df = pd.read_excel('data/Copy of reversible lambda for Jihoon.xlsx', comment='#', )
df
Out[5]:
drug condition name days after starting drug treatment l1 l2
0 drug off Control 0 47.372417 -32.483135
1 drug on D3 3 68.983442 36.353281
2 drug on D8 8 49.936484 39.770341
3 drug on D13 13 31.444846 38.371612
4 drug on D21 21 -8.357498 31.326270
5 drug on D29 29 -35.907250 19.199379
6 drug on D33 33 -41.172003 21.773860
7 drug on D38 38 -47.742299 10.481910
8 drug on D59 59 -41.966089 18.962689
9 drug on D29 29 -35.907250 19.199379
10 drug off dr29.D4 33 -45.236981 -12.113332
11 drug off dr29.D10 39 -35.883245 -19.890666
12 drug off dr29.D15 44 -12.973426 -32.785915
13 drug off dr29.D17 46 -2.289392 -37.153589
14 drug off dr29.D30 59 36.139257 -38.247608
15 drug off dr29.D35 64 38.022285 -37.018481
In [6]:
df['l1'] = df['l1'] + 50.0
df['l2'] = df['l2'] + 50.0
In [7]:
const_drug = df[:9]
pulse_drug = df.drop(df.index[6:10])
pulse_drug = pulse_drug.drop(pulse_drug.index[0:5])
const_drug = const_drug.sort('days after starting drug treatment')
pulse_drug = pulse_drug.sort('days after starting drug treatment')
const_drug_alt = const_drug
const_drug_alt
C:\Users\Jihoon\Anaconda3\lib\site-packages\ipykernel\__main__.py:4: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)
C:\Users\Jihoon\Anaconda3\lib\site-packages\ipykernel\__main__.py:5: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)
Out[7]:
drug condition name days after starting drug treatment l1 l2
0 drug off Control 0 97.372417 17.516865
1 drug on D3 3 118.983442 86.353281
2 drug on D8 8 99.936484 89.770341
3 drug on D13 13 81.444846 88.371612
4 drug on D21 21 41.642502 81.326270
5 drug on D29 29 14.092750 69.199379
6 drug on D33 33 8.827997 71.773860
7 drug on D38 38 2.257701 60.481910
8 drug on D59 59 8.033911 68.962689
In [39]:
# Add a few dummy data points at the end to try to induce a steady state equilibrium and deter oscillation
# to pulsed case
pulse_drug = pulse_drug.append(pd.DataFrame([['drug off', 'dr29.D40', 69, 88.022285, 12.981519], \
                                            ['drug off', 'dr29.D45', 74, 88.022285, 12.981519], \
                                            ['drug off', 'dr29.D50', 79, 88.022285, 12.981519]], \
                               columns=['drug condition', 'name', 'days after starting drug treatment', 'l1', 'l2']))
pulse_drug = pulse_drug.reset_index()
pulse_drug = pulse_drug.drop('index', 1)
const_drug = pulse_drug
const_drug
Out[39]:
drug condition name days after starting drug treatment l1 l2
0 drug on D29 29 14.092750 69.199379
1 drug off dr29.D4 33 4.763019 37.886668
2 drug off dr29.D10 39 14.116755 30.109334
3 drug off dr29.D15 44 37.026574 17.214085
4 drug off dr29.D17 46 47.710608 12.846411
5 drug off dr29.D30 59 86.139257 11.752392
6 drug off dr29.D35 64 88.022285 12.981519
7 drug off dr29.D40 69 88.022285 12.981519
8 drug off dr29.D45 74 88.022285 12.981519
9 drug off dr29.D50 79 88.022285 12.981519
In [40]:
# Theoretical model for populations
def popDynamics(p, d):
    data_times, (c01, c02) = d
    data_times = np.array(data_times)
    a1, a2, a3, alpha1, beta1, alpha2, s = p
    euler_step = 0.5
    td = np.arange(data_times[0], data_times[-1] + euler_step, euler_step)
    l_theor_full = np.zeros([len(td), len(d[1])]) 
    l_theor_full[0] = np.array([[c01, c02]])
    l_theor = np.array([[c01, c02]])
    
    # Temporarily? changed the diff eqs such that a1 and a2 are the same, in order to make the resulting vector field
    # conservative --> more biologically plausible
    for i in range(1, len(td)):
        l_theor_full[i][0] = l_theor_full[i-1][0] + euler_step * (alpha1 - beta1 * l_theor_full[i-1][0] \
                                                                 - a1 * l_theor_full[i-1][1])
        l_theor_full[i][1] = l_theor_full[i-1][1] + euler_step * (alpha2 + a2 * l_theor_full[i-1][0] \
                                                                 - a3 * l_theor_full[i-1][1])
        if td[i] in data_times:
            l_theor = np.append(l_theor, [[l_theor_full[i][0], l_theor_full[i][1]]], axis=0)
    return l_theor_full, l_theor
    
def log_posterior(p, d, l):
    '''
    Assuming uniform priors for all of the parameters and Jeffreys prior for sigma. Assuming Gaussian 
    distributed measurements in the experiment. 
    '''
    data_times, (c01, c02) = d
    a1, a2, a3, alpha1, beta1, alpha2, s = p
    # Zero probability of having non-positive value of stdev
    if s <= 0:
        return -np.inf
    _, pops_theor = popDynamics(p, d)
    # Zero probability of having absurd population values
    for i in range(len(pops_theor)):
        for j in range(2):
            if np.abs(pops_theor[i][j]) > 300.0 or np.isnan(pops_theor[i][j]) or np.abs(pops_theor[i][j]) < 0:
                return -np.inf
    for val in p:
        if val < 0:
            return -np.inf
    #if alpha1 < 0 or alpha2 < 0 or a2 < 0 or a1 < 0:
        #return -np.inf
    if a1 > 1 or a2 > 1 or a3 > 1:
        return -np.inf

    # Nyquist frequency: 1 day^-1 so keep frequency below or equal to that
    if 1 / (2 * np.pi) * np.sqrt(np.sqrt(a1**2.0 * a3**2.0) + np.sqrt(beta1**2.0 * a2**2.0)) > 1.0:
        return -np.inf
    return -(2 * len(data_times) + 1) * np.log(s) - 0.5 / s**2.0 * (np.sum((l[:,1]-pops_theor[:,1])**2.0) \
                                             + np.sum((l[:,0]-pops_theor[:,0])**2.0))
In [43]:
n_dim = 7 # 7 parameters in the model
n_walkers = 1000 # Number of MCMC walkers
n_burn = 500 # 500 burn-in steps
n_steps = 5000 # Total number of steps after burn-in

# Seeding random number generator 
np.random.seed(1)
In [44]:
# p0[i, j] is the starting point for the ith walk for the jth variable
p0 = np.empty((n_walkers, n_dim))

for i in range(n_dim - 1):
    p0[:, i] = np.random.uniform(0, 1, n_walkers)
p0[:, n_dim-1] = np.random.exponential(0.1, n_walkers) #sigma
In [45]:
# Creating a MCMC sampler
emp_data = np.array([const_drug['l1'], const_drug['l2']])
emp_data = np.transpose(emp_data)
sampler = emcee.EnsembleSampler(n_walkers, n_dim, log_posterior, args=((const_drug['days after starting drug treatment'], \
        emp_data[0]), emp_data,), threads=1) 
        # Using 1 core since multithreading seems to be broken
In [10]:
# Do burn-in
pos, prob, state = sampler.run_mcmc(p0, n_burn, storechain=False)

#Actual run
for i, result in enumerate(sampler.sample(pos, lnprob0=prob, iterations=n_steps)):
    if (i+1) % 100 == 0:
        print("{0:5.1%}".format(float(i) / n_steps))
 2.0%
 4.0%
 6.0%
 8.0%
10.0%
12.0%
14.0%
16.0%
18.0%
20.0%
22.0%
24.0%
26.0%
28.0%
30.0%
32.0%
34.0%
36.0%
38.0%
40.0%
42.0%
44.0%
46.0%
48.0%
50.0%
52.0%
54.0%
56.0%
58.0%
60.0%
62.0%
64.0%
66.0%
68.0%
70.0%
72.0%
74.0%
76.0%
78.0%
80.0%
82.0%
84.0%
86.0%
88.0%
90.0%
92.0%
94.0%
96.0%
98.0%
100.0%
In [11]:
fig, ax = plt.subplots(7, 1, sharex=True)
for i in range(6):
    ax[i].plot(sampler.chain[0,:,i], 'k-', lw=0.4)
    ax[i].plot([0, n_steps-1], 
             [sampler.chain[0,:,i].mean(), sampler.chain[0,:,i].mean()], 'r-', lw=0.4)

ax[6].set_xlabel('sample number')
Out[11]:
<matplotlib.text.Text at 0x22e8e85fac8>
In [12]:
# Get the index of the most probable parameter set
max_ind = np.argmax(sampler.flatlnprobability)

# Pull out values.
a1, a2, a3, alpha1, beta1, alpha2, s = sampler.flatchain[max_ind,:]
# Calculate errors
a1_e, a2_e, a3_e, alpha1_e, beta1_e, alpha2_e, s_e = sampler.flatchain.std(axis=0)

# Print the results
print('Most probable parameter values:' + \
      '\na1: ' + str(a1) + ' +- ' + str(a1_e) + \
      '\na2: ' + str(a2) + ' +- ' + str(a2_e) + \
      '\na3: ' + str(a3) + ' +- ' + str(a3_e) + \
      '\nalpha1: ' + str(alpha1) + ' +- ' + str(alpha1_e) + \
      '\nbeta1: ' + str(beta1) + ' +- ' + str(beta1_e) + \
      '\nalpha2: ' + str(alpha2) + ' +- ' + str(alpha2_e) + \
      '\ns: ' + str(s) + ' +- ' + str(s_e))
Most probable parameter values:
a1: 0.247519913164 +- 0.237015464185
a2: 0.0160078926305 +- 0.24349486565
a3: 0.109295606836 +- 0.314391329252
alpha1: 11.6654030819 +- 216.335085565
beta1: 0.090416552161 +- 1.94643166363
alpha2: 0.123918029102 +- 6.97259287742
s: 2.68210259427 +- 6.22285059874
In [261]:
_ = plt.plot(td, popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, [0.0, 0.0]))[0], alpha=0.8)
In [13]:
_ = plt.scatter(const_drug['days after starting drug treatment'], const_drug['l2'], color='green')
_ = plt.scatter(const_drug['days after starting drug treatment'], const_drug['l1'], color='blue')
euler_step = 0.5
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
               np.array(const_drug['days after starting drug treatment'])[-1] + euler_step, euler_step)

plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(td, popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, emp_data[0]))[0], lw=0.6, alpha=0.8)
_ = plt.plot(const_drug['days after starting drug treatment'], \
             popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (const_drug['days after starting drug treatment'], emp_data[0]))[1], lw=1.5)
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
C:\Users\Jihoon\Anaconda3\lib\site-packages\matplotlib\cbook.py:137: MatplotlibDeprecationWarning: The set_color_cycle attribute was deprecated in version 1.5. Use set_prop_cycle instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [14]:
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
               np.array(const_drug['days after starting drug treatment'])[-1] + 100 * euler_step, euler_step)
_ = plt.plot(td, popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, emp_data[0]))[0])
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
np.savetxt('data/pulse_drug_results.csv', popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, emp_data[0]))[0], \
          delimiter=',')
In [15]:
# Plotting the probability distributions of the most probable parameter values
#fig = corner.corner(sampler.flatchain, labels=[r'$a_1$', r'$a_2$', r'$a_3$', r'$\alpha_1$', r'$\beta_1$', \
#                                                  r'$\alpha_2$', r'$\sigma$'], bins=100, lw=2)
#fig.savefig("data/corner_pulse.png")

Now we do the same for the constant drug (no release):

In [52]:
euler_step = 0.5
emp_data_alt = np.array([const_drug_alt['l1'], const_drug_alt['l2']])
emp_data_alt = np.transpose(emp_data_alt)

sampler_const = emcee.EnsembleSampler(n_walkers, n_dim, log_posterior, args=((const_drug_alt['days after starting drug treatment'], \
        emp_data_alt[0]), emp_data_alt,), threads=1) 
        # Using 1 core since multithreading seems to be broken
In [418]:
# Do burn-in
pos, prob, state = sampler_const.run_mcmc(p0, n_burn, storechain=False)

for i, result in enumerate(sampler_const.sample(pos, lnprob0=prob, iterations=n_steps)):
    if (i+1) % 100 == 0:
        print("{0:5.1%}".format(float(i) / n_steps))

# Get the index of the most probable parameter set
max_ind = np.argmax(sampler_const.flatlnprobability)

# Pull out values.
a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc = sampler_const.flatchain[max_ind,:]
# Calculate errors
a1c_e, a2c_e, a3c_e, alpha1c_e, beta1c_e, alpha2c_e, sc_e = sampler_const.flatchain.std(axis=0)

# Print the results
print('Most probable parameter values:' + \
      '\na1: ' + str(a1c) + ' +- ' + str(a1c_e) + \
      '\na2: ' + str(a2c) + ' +- ' + str(a2c_e) + \
      '\na3: ' + str(a3c) + ' +- ' + str(a3c_e) + \
      '\nalpha1: ' + str(alpha1c) + ' +- ' + str(alpha1c_e) + \
      '\nbeta1: ' + str(beta1c) + ' +- ' + str(beta1c_e) + \
      '\nalpha2: ' + str(alpha2c) + ' +- ' + str(alpha2c_e) + \
      '\ns: ' + str(sc) + ' +- ' + str(sc_e))

td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
               np.array(const_drug_alt['days after starting drug treatment'])[-1] + 100 * euler_step, euler_step)
_ = plt.plot(td, popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (td, emp_data_alt[0]))[0])
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
np.savetxt('data/const_drug_results.csv', popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (td, emp_data_alt[0]))[0], \
          delimiter=',')
 2.0%
 4.0%
 6.0%
 8.0%
10.0%
12.0%
14.0%
16.0%
18.0%
20.0%
22.0%
24.0%
26.0%
28.0%
30.0%
32.0%
34.0%
36.0%
38.0%
40.0%
42.0%
44.0%
46.0%
48.0%
50.0%
52.0%
54.0%
56.0%
58.0%
60.0%
62.0%
64.0%
66.0%
68.0%
70.0%
72.0%
74.0%
76.0%
78.0%
80.0%
82.0%
84.0%
86.0%
88.0%
90.0%
92.0%
94.0%
96.0%
98.0%
100.0%
Most probable parameter values:
a1: 0.338521518145 +- 0.261579331293
a2: 0.0759537496773 +- 0.28456167881
a3: 0.399963622611 +- 0.20527522408
alpha1: 24.1813265259 +- 47.5655886632
beta1: 0.000262834618533 +- 1.72851655608
alpha2: 28.6176624135 +- 18.5808178781
s: 6.06207055234 +- 6.74430878985
In [419]:
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l2'], color='green')
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l1'], color='blue')
td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
               np.array(const_drug_alt['days after starting drug treatment'])[-1] + euler_step, euler_step)
plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(td, popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (td, emp_data_alt[0]))[0], lw=0.6, alpha=0.8)
_ = plt.plot(const_drug_alt['days after starting drug treatment'], \
             popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (const_drug_alt['days after starting drug treatment'], emp_data_alt[0]))[1], lw=1.5)
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
C:\Users\Jihoon\Anaconda3\lib\site-packages\matplotlib\cbook.py:137: MatplotlibDeprecationWarning: The set_color_cycle attribute was deprecated in version 1.5. Use set_prop_cycle instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [18]:
# Plotting the probability distributions of the most probable parameter values CONSTANT DRUG
#fig = corner.corner(sampler_const.flatchain, labels=[r'$a_1$', r'$a_3$', r'$\alpha_1$', r'$\beta_1$', \
#                                                  r'$\alpha_2$', r'$\sigma$'], bins=100, lw=2)
#fig.savefig('data/corner_const.png')
In [88]:
def distances(p, l1Domain, l2Domain, td, pulseTrue):
    # Calculating trajectory length from start to steady state
    # l1Domain and l2Domain are the starting domains of l1 and l2. Must be equally sized
    # Overlays the trajectory using empirical starting points as well, with z-value equal to distance left to cover.
    if pulseTrue:
        dist_emp_data = emp_data[0]
    else:
        dist_emp_data = emp_data_alt[0]
    a1, a2, a3, alpha1, beta1, alpha2, s = p
    euler_step = 0.5
    dist_comb = 0
    l_theor_full = np.zeros([len(td), 2, len(l1Domain), len(l2Domain)]) 
    l_theor_full[0] = np.meshgrid(l1Domain, l2Domain)
    for i in range(1, len(td)):
        l_theor_full[i][0] = l_theor_full[i-1][0] + euler_step * (alpha1 - beta1 * l_theor_full[i-1][0] \
                                                                 - a1 * l_theor_full[i-1][1])
        l_theor_full[i][1] = l_theor_full[i-1][1] + euler_step * (alpha2 + a2 * l_theor_full[i-1][0] \
                                                                 - a3 * l_theor_full[i-1][1])
        dist_1 = l_theor_full[i][0] - l_theor_full[i-1][0]
        dist_2 = l_theor_full[i][1] - l_theor_full[i-1][1]
        dist_comb = dist_comb + np.sqrt(dist_1**2.0 + dist_2**2.0)
    
    l1Domain, l2Domain = np.meshgrid(l1Domain, l2Domain)
    pops = popDynamics(p, (td, dist_emp_data))[0]
    dist_traj_indiv = np.sqrt(np.diff([val[0] for val in pops])**2.0 + np.diff([val[1] for val in pops])**2.0)
    dist_traj_indiv = np.append(dist_traj_indiv, dist_traj_indiv[-1])
    dist_traj = np.zeros(len(dist_traj_indiv))
    for i in range(len(dist_traj_indiv)):
        dist_traj[i] = np.sum(dist_traj_indiv) - np.sum(dist_traj_indiv[:i])
    #dist_traj = dist_traj[::-1]
    
    UN = alpha1 - beta1 * l1Domain - a1 * l2Domain
    VN = alpha2 + a2 * l1Domain - a3 * l2Domain
    fig = plt.figure(figsize=(12,9))
    ax = fig.gca(projection='3d')
    pulse_surf = ax.plot_surface(l1Domain, l2Domain, dist_comb, cstride=10, rstride=10, cmap=cm.jet, \
                                 linewidth=0, antialiased=False, alpha=0.6)
    fig.colorbar(pulse_surf, shrink=0.5, aspect=5)
    pulse_surf = ax.quiver(l1Domain[::20, ::20], l2Domain[::20, ::20], dist_comb[::20, ::20], UN[::20, ::20], \
                           VN[::20, ::20], 0, length=1.0, linewidths=1.5, arrow_length_ratio=0.6, color='black', pivot='middle')
    pulse_surf = ax.plot([val[0] for val in pops], [val[1] for val in pops], dist_traj, lw=1.4, color='white', \
                        label='From exp. starting point')
    pulse_surf = ax.set_xlabel(r'$\lambda_1$ starting value')
    pulse_surf = ax.set_ylabel(r'$\lambda_2$ starting value')
    pulse_surf = ax.set_zlabel('Total distance of trajectory to steady state')
    if pulseTrue == False:
        ax.view_init(azim=200, elev=30)
        ax.set_title('Constant drug condition trajectory lengths')
    else:
        ax.set_title('Pulse drug condition trajectory lengths')
        ax.view_init(azim=-40, elev=30)
    return np.min(np.min(dist_comb))
In [82]:
s = 2.0
euler_step = 0.5
# Pulse drug
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
               np.array(const_drug['days after starting drug treatment'])[-1] + euler_step, euler_step)
l1Domain = np.arange(-30, 210, 1)
l2Domain = np.arange(-30, 210, 1)
distances((a1, a2, a3, alpha1, beta1, alpha2, s), l1Domain, l2Domain, td, True)
Out[82]:
0.79511751619249482
In [89]:
sc = 2.0
# Constant drug
td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
               np.array(const_drug_alt['days after starting drug treatment'])[-1] + euler_step, euler_step)
distances((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), l1Domain, l2Domain, td, False)
Out[89]:
0.75868886766626908
In [407]:
def distances_2d(p, l1Domain, l2Domain, td, pulseTrue):
    # Calculating trajectory length from start to steady state
    # l1Domain and l2Domain are the starting domains of l1 and l2. Must be equally sized
    # Overlays the trajectory using empirical starting points as well
    if pulseTrue:
        dist_emp_data = emp_data[0]
    else:
        dist_emp_data = emp_data_alt[0]
    a1, a2, a3, alpha1, beta1, alpha2, s = p
    euler_step = 0.5
    dist_comb = 0
    l_theor_full = np.zeros([len(td), 2, len(l1Domain), len(l2Domain)]) 
    l_theor_full[0] = np.meshgrid(l1Domain, l2Domain)
    for i in range(1, len(td)):
        l_theor_full[i][0] = l_theor_full[i-1][0] + euler_step * (alpha1 - beta1 * l_theor_full[i-1][0] \
                                                                 - a1 * l_theor_full[i-1][1])
        l_theor_full[i][1] = l_theor_full[i-1][1] + euler_step * (alpha2 + a2 * l_theor_full[i-1][0] \
                                                                 - a3 * l_theor_full[i-1][1])
        dist_1 = l_theor_full[i][0] - l_theor_full[i-1][0]
        dist_2 = l_theor_full[i][1] - l_theor_full[i-1][1]
        dist_comb = dist_comb + np.sqrt(dist_1**2.0 + dist_2**2.0)
    
    l1Domain, l2Domain = np.meshgrid(l1Domain, l2Domain)
    pops = popDynamics(p, (td, dist_emp_data))[0]
    dist_traj_indiv = np.sqrt(np.diff([val[0] for val in pops])**2.0 + np.diff([val[1] for val in pops])**2.0)
    dist_traj_indiv = np.append(dist_traj_indiv, dist_traj_indiv[-1])
    dist_traj = np.zeros(len(dist_traj_indiv))
    for i in range(len(dist_traj_indiv)):
        dist_traj[i] = np.sum(dist_traj_indiv) - np.sum(dist_traj_indiv[:i])
    
    UN = alpha1 - beta1 * l1Domain - a1 * l2Domain
    VN = alpha2 + a2 * l1Domain - a3 * l2Domain
    fig = plt.figure(figsize=(12,9))
    ax = fig.gca()
    pulse_surf = plt.plot([val[0] for val in pops], [val[1] for val in pops], lw=1.4, color='white')
    pulse_surf = plt.pcolormesh(l1Domain, l2Domain, dist_comb, cmap=cm.jet, \
                                 linewidth=0, antialiased=False, alpha=0.6)
    fig.colorbar(pulse_surf, shrink=0.5, aspect=5, label='Trajectory distance')
    pulse_surf = plt.quiver(l1Domain[::20, ::20], l2Domain[::20, ::20], UN[::20, ::20], \
                           VN[::20, ::20], color='black', pivot='middle')
    pulse_surf = plt.xlabel(r'$\lambda_1$ starting value')
    pulse_surf = plt.ylabel(r'$\lambda_2$ starting value')
    if pulseTrue == False:
        plt.title('Constant drug condition trajectory lengths')
    else:
        plt.title('Pulse drug condition trajectory lengths')
    return np.min(np.min(dist_comb))
In [408]:
# Pulse drug
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
               np.array(const_drug['days after starting drug treatment'])[-1] + euler_step, euler_step)
l1Domain = np.arange(-30, 210, 1)
l2Domain = np.arange(-30, 210, 1)
distances_2d((a1, a2, a3, alpha1, beta1, alpha2, s), l1Domain, l2Domain, td, True)
Out[408]:
0.64723036584159976
In [409]:
# Constant drug
# Pulse drug
td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
               np.array(const_drug_alt['days after starting drug treatment'])[-1] + euler_step, euler_step)
l1Domain = np.arange(-30, 210, 1)
l2Domain = np.arange(-30, 210, 1)
distances_2d((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), l1Domain, l2Domain, td, False)
Out[409]:
0.21998704127846155
In [27]:
# More analysis of previously found model fitting
df_pulse = pd.read_csv('data/pulse_drug_results.csv', comment='#', )
df_const = pd.read_csv('data/const_drug_results.csv', comment='#',)
df_pulse.head()
Out[27]:
time (cycles) l1 l2
0 29.0 14.10 69.2
1 29.5 10.70 65.6
2 30.0 7.95 62.2
3 30.5 5.74 58.9
4 31.0 4.02 55.8
In [129]:
_ = plt.scatter(pulse_drug['days after starting drug treatment'], pulse_drug['l2']-50.0, color='green')
_ = plt.scatter(pulse_drug['days after starting drug treatment'], pulse_drug['l1']-50.0, color='blue')
plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(df_pulse['time (cycles)'], df_pulse['l1']-50.0, df_pulse['time (cycles)'], df_pulse['l2']-50.0)
_ = plt.legend([r'$\lambda_1$', r'$\lambda_2$'], loc='upper right')
_ = plt.xlabel('Day')
_ = plt.ylabel('population')
_ = plt.show()
C:\Users\Jihoon\Anaconda3\lib\site-packages\matplotlib\cbook.py:137: MatplotlibDeprecationWarning: The set_color_cycle attribute was deprecated in version 1.5. Use set_prop_cycle instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [130]:
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l2']-50.0, color='green')
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l1']-50.0, color='blue')
plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(df_const['time (cycle)'], df_const['l1']-50.0, df_const['time (cycle)'], df_const['l2']-50.0,)
_ = plt.legend([r'$\lambda_1$', r'$\lambda_2$'], loc='upper right')
_ = plt.xlabel('Day')
_ = plt.ylabel('population')
_ = plt.show()
C:\Users\Jihoon\Anaconda3\lib\site-packages\matplotlib\cbook.py:137: MatplotlibDeprecationWarning: The set_color_cycle attribute was deprecated in version 1.5. Use set_prop_cycle instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)
In [9]:
df_landscape = pd.read_excel('data/input_reversible_landscape.xlsx', comment='#')
df_landscape.head()
Out[9]:
xi yi z
0 -70 -70 171510.264029
1 -70 -63 163376.703844
2 -70 -56 155090.974248
3 -70 -49 146502.848963
4 -70 -42 137614.225620
In [112]:
xi, yi = np.meshgrid(df_landscape['xi'].unique() + 50.0, df_landscape['yi'].unique() + 50.0)
z = np.zeros((len(xi), len(yi)))
for i in range(len(xi)):
    for j in range(len(yi)):
        z[i][j] = df_landscape['z'][j*len(xi) + i] / 1000.0

a1 = 0.248
a2 = 0.016
a3 = 0.109
alpha1 = 11.665
beta1 = 0.090
alpha2 = 0.124
UN = alpha1 - beta1 * xi - a1 * yi
VN = alpha2 + a2 * xi - a3 * yi
        
# Plotting landscape
fig = plt.figure(figsize=(12, 8))
ax = fig.gca(projection='3d')
# Create light source object.
ls = LightSource(azdeg=285, altdeg=70, hsv_min_val=0.6, hsv_min_sat=0.8)
# Shade data, creating an rgb array.
rgb = ls.shade(z, plt.cm.rainbow)
pulse_surf = ax.plot_surface(xi-50, yi-50, z, cstride=1, rstride=1, facecolors=rgb, \
                             linewidth=0, antialiased=True, alpha=1, shade=True)

pulse_surf = ax.quiver(xi-50, yi-50, z, UN, \
                           VN, 0, length=0.65, linewidths=1.5, arrow_length_ratio=0.6, color='black', pivot='middle')
pulse_surf = ax.set_xlabel(r'$\lambda_1$ starting value')
pulse_surf = ax.set_ylabel(r'$\lambda_2$ starting value')
pulse_surf = ax.set_zlabel(r'Free energy ($\times10^3$)')
ax.view_init(azim=240, elev=60)
#ax.dist=0.1
plt.title('Pulse drug condition energy landscape')
Out[112]:
<matplotlib.text.Text at 0x25408646f60>
In [113]:
xi, yi = np.meshgrid(df_landscape['xi'].unique() + 50.0, df_landscape['yi'].unique() + 50.0)
z = np.zeros((len(xi), len(yi)))
for i in range(len(xi)):
    for j in range(len(yi)):
        z[i][j] = df_landscape['z'][j*len(xi) + i] / 1000.0

a1c = 0.339
a2c = 0.076
a3c = 0.400
alpha1c = 24.181
beta1c = 0.0003
alpha2c = 28.618
UN = alpha1c - beta1c * xi - a1c * yi
VN = alpha2c + a2c * xi - a3c * yi
        
# Plotting landscape
fig = plt.figure(figsize=(12, 8))
ax = fig.gca(projection='3d')
# Create light source object.
ls = LightSource(azdeg=285, altdeg=70, hsv_min_val=0.6, hsv_min_sat=0.8)
# Shade data, creating an rgb array.
rgb = ls.shade(z, plt.cm.rainbow)
pulse_surf = ax.plot_surface(xi-50, yi-50, z, cstride=1, rstride=1, facecolors=rgb, \
                             linewidth=0, antialiased=True, alpha=1, shade=True)

pulse_surf = ax.quiver(xi-50, yi-50, z, UN, \
                           VN, 0, length=0.65, linewidths=1.5, arrow_length_ratio=0.6, color='black', pivot='middle')
pulse_surf = ax.set_xlabel(r'$\lambda_1$ starting value')
pulse_surf = ax.set_ylabel(r'$\lambda_2$ starting value')
pulse_surf = ax.set_zlabel(r'Free energy ($\times10^3$)')
ax.view_init(azim=240, elev=60)
#ax.dist=0.1
plt.title('Constant drug condition energy landscape')
Out[113]:
<matplotlib.text.Text at 0x2540c3d2a90>
In [ ]: